
#include "interpolation.h"

// 基类析构函数必须有实现
PolynomialInterpolation::~PolynomialInterpolation() {}

// 返回插值多项式在x的值
double NewtonInterpolation::evaluate(double x) const
{
    double result = 0.0, pik = 1.0;
    for (int i = 0; i <= n; i++)
    {
        result += divided_differences_table[i][i] * pik;
        pik *= (x - interpolation_points[i]);
    }
    return result;
}

// 计算差商表
void NewtonInterpolation::calculate()
{
    for (int j = 1; j <= n; j++)
    {
        for (int i = j; i <= n; i++)
        {
            divided_differences_table[i][j] =
                (divided_differences_table[i][j - 1] - divided_differences_table[i - 1][j - 1]) / (interpolation_points[i] - interpolation_points[i - j]);
        }
    }
    return;
}

// 打印区间[left，right]内的插值函数
void NewtonInterpolation::print(const double left, const double right, const double grid) const
{
    for (double x = left; x <= right; x += grid)
    {
        std::cout << x << " " << evaluate(x) << std::endl;
    }
    return;
}

// 将区间[left，right]内的插值函数输出到txt文件中
void NewtonInterpolation::output_file(const char *filename, const double left, const double right,
                                      const double grid) const
{
    std::ofstream outputFile(filename); // 打开输出文件

    std::streambuf *originalCout = std::cout.rdbuf(); // 保存原始的 std::cout 缓冲区

    std::cout.rdbuf(outputFile.rdbuf()); // 重定向 std::cout 到输出文件

    print(left, right, grid);

    std::cout.rdbuf(originalCout); // 恢复原始的 std::cout 缓冲区

    outputFile.close(); // 关闭文件
}

NewtonInterpolation::NewtonInterpolation(const std::vector<double> &ips, int np, double (*f)(double))
    : PolynomialInterpolation{ips, np},
      divided_differences_table{std::vector<std::vector<double>>(np + 1, std::vector<double>(np + 1, 0.0))}
{
    // 预处理差商表
    for (int i = 0; i <= n; i++)
    {
        divided_differences_table[i][0] = f(ips[i]);
    }
    calculate();
}

NewtonInterpolation::NewtonInterpolation(const std::vector<double> &ips, int np, const std::vector<double> &ivs)
    : PolynomialInterpolation{ips, np},
      divided_differences_table{std::vector<std::vector<double>>(np + 1, std::vector<double>(np + 1, 0.0))}
{
    // 预处理差商表
    for (int i = 0; i <= n; i++)
    {
        divided_differences_table[i][0] = ivs[i];
    }
    calculate();
}

ChebyshevInterpolation::ChebyshevInterpolation(int np, double (*f)(double))
{
    n = np;
    interpolation_points = std::vector<double>(n + 1);
    divided_differences_table = std::vector<std::vector<double>>(n + 1, std::vector<double>(n + 1, 0.0));
    for (int i = 0; i <= n; i++)
    {
        interpolation_points[i] = cos((i + 0.5) / (n + 1) * M_PI);    // 生成插值点
        divided_differences_table[i][0] = f(interpolation_points[i]); // 预处理差商表
    }
    calculate(); // 按照牛顿插值计算差商表
}

double HermiteInterpolation::diff(double x, double eps) const
{
    return (evaluate(x + eps) - evaluate(x - eps)) / (2 * eps);
}

// 计算差商表
void HermiteInterpolation::calculate()
{
    for (int j = 1; j <= n; j++)
    {
        for (int i = j; i <= n; i++)
        {
            // 若分母为0,则此处为重结点差商，已存在导数信息，无需计算
            if (interpolation_points[i] - interpolation_points[i - j] == 0)
                continue;
            divided_differences_table[i][j] =
                (divided_differences_table[i][j - 1] - divided_differences_table[i - 1][j - 1]) /
                (interpolation_points[i] - interpolation_points[i - j]);
        }
    }
    return;
}

HermiteInterpolation::HermiteInterpolation(const std::vector<double> &ips, const std::vector<int> &order,
                                           const std::vector<std::vector<double>> &ivs)
{
    int s = order.size() - 1;
    n = accumulate(order.begin(), order.end(), s);
    interpolation_points = std::vector<double>(n + 1);
    divided_differences_table = std::vector<std::vector<double>>(n + 1, std::vector<double>(n + 1, 0.0));
    int t = 0;
    // 预处理差商表，并处理插值点
    for (int i = 0; i <= s; i++)
    {
        for (int j = 0; j <= order[i]; j++)
        {
            int factorial;
            interpolation_points[t] = ips[i];
            for (int k = 0; k <= j; k++)
            {
                if (k == 0)
                    factorial = 1;
                else
                    factorial *= k;
                divided_differences_table[t][k] = ivs[i][k] / factorial;
            }
            t++;
        }
    }
    calculate();
}
